% 快速傅里叶变换的演示文档
% * 查看信号的周期
% * 快速傅里叶变换函数
% 自定义函数：
%	* myFFT：	快速傅里叶变换函数
%	* myFFTshow：	带绘图分析的快速傅里叶变换函数	
% 作者：sid（cnqdcyq@163.com）
% 时间：2021/7/28
% 参考：[快速傅里叶变换-MATHWORKS](https://ww2.mathworks.cn/help/matlab/ref/fft.html)
t = transpose(1:100);


% 定义信号
Fs = 1000; % 取样频率,单位时间取样的次数
T = 1/Fs; % 取样周期
L = 1500; % 信号长度
t = transpose( (0:L-1) * T ); % 时间向量
Signal = 0.7 * sin(2*pi*50*t) + sin(2*pi*120*t);
Noise =  2*randn(size(t));
Input = Signal + Noise;
plot(t(1:50),Input(1:50))


% 进行信号分解
[f,P1] = myFFT(Input,Fs);


% 展示周期信息
myFFTshow(Input,Fs);

% 自定义函数区
function [f,P1] = myFFT(input, Fs)
	% 快速傅里叶变换，输入时间序列，分解为频域
	% 输入变量：
	% 	input：待处理的信号
	%	Fs：	取样频率
	L = length(input);
	Y = fft(input);
	P2 = abs(Y/L); % 计算双侧频谱
	P1 = P2(1:L/2+1); 
	P1(2:end-1) = 2*P1(2:end-1); % 计算单侧频谱

	f = Fs * (0:(L/2)) / L;
end

function [f,P1] = myFFTshow(input,Fs)
	% 快速傅里叶变换，输入时域信号转换为频域信号
	% 输入变量：
	% 	input：	待处理的信号，列向量
	%	Fs：	取样频率
	% 输出变量：
	%	f：	频率
	%	P1:	振幅
	[f,P1] = myFFT(input,Fs);
	plot(f,P1)
	title("单侧振幅频谱图")
	xlabel("f(Hz)");
	ylabel("|P1(f)|");

end
